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two-dimensional couplings 

Pinaki Sengupta,^'Q Weihong Zheng, and Rajiv R. P. Singh^'0 

'Department of Physics, University of California, Davis, California 95616 
"^School of Physics, University of New South Wales, Sydney, NSW 2052, Australia 

(Dated: February 2, 2008) 

Motivated by the geometry of spins in the material CaCu203, we study a two-layer, spin-half 
Heisenberg model, with nearest-neighbor exchange couplings J and aJ along the two axes in the 
plane and a coupling J± perpendicular to the planes. We study these class of models using the 
Stochastic Series Expansion (SSE) Quantum Monte Carlo simulations at finite temperatures and 
series expansion methods at T = 0. The critical value of the interlayer coupling, J^, separating the 
Neel ordered and disordered ground states, is found to follow very closely a square root dependence 
on a. Both T = and finite-temperature properties of the model are presented. 

PACS numbers: PACS: 75.40.Gb, 75.40.Mg, 75.10.Jm, 75.30.Ds 
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I. INTRODUCTION 

In recent years, the geometrical arrangement of atoms 
carrying spin and those mediating interactions between 
them, in newly discovered magnetic materials, has been 
a dominant theme in the field of quantum magnetism, 
and has inspired many interesting models and theoriesii. 
In a class of cuprate materials, the arrangement of cop- 
per and oxygen atoms in the copper-oxide planes leads 
to nearly ideal realizations of ID spin-chains, of spin- 
ladders as well as of simple and decorated square-lattice 
systems^. The variations in geometry, in other recently 
discovered materials, have led to several magnetically or- 
dered, spin-gap and dimerized phases'*. These materials 
have allowed the exploration of dimensional crossovers, 
unusual excitations and quantum critical phenomenal 
through detailed comparisons between experiments and 
theory. 

Many classes of materials involve weakly coupled spin- 
half chains. From a theoretical point of view, these are 
ideal systems for finding and testing novel phenomena as 
our understanding of the spin-half chain is rather com- 
plete thanks to Bcthe ansatz methods, quantum field- 
theory methods and numerical techniques^'^'^. However, 
different ways of coupling the chains together can lead to 
widely different behaviors. If two Heisenberg chains are 
coupled weakly together, it is known to lead to a disor- 
dered ground state and a gap in the spin-excitation spec- 
tra. These systems have been called spin-ladders and are 
considered generic examples of spin-gap phenomena^*iS. 
On the other hand, if an array of spin-chains, say ar- 
ranged parallel to each other in a two-dimensional plane, 
are weakly coupled together in an unfrustrated manner, 
they develop conventional Neel order—. If the inter- 
chain couplings are strongly frustrated, one can even find 
a carryover of exotic one-dimensional physics to higher 
dimensionsiSii^. 

In this paper, we are motivated by the material 
CaCu203, which has a geometry very similar to the spin- 



temperaturesii. The largest exchange coupling in these 
materials is of order 2000/ir and leads to spin-half chains. 
These chains are weakly coupled in two-different ways. 
Along one axis perpendicular to the chains, one has ef- 
fectively a two-leg ladder, which are then nearly decou- 
pled from the rest of the system, due to the geometry 
of copper-oxygen bond angles. Along the other axis per- 
pendicular to the chains, the coupling between the chains 
leads to a square-lattice of spins-^. Thus the material 
has competing tendencies to Neel order and to develop a 
gap in the spin excitation spectra. Neel order is foundi^ 
in bulk materials with T/v = 25K. In a Science paper, 
Schon et. al— had argued that thin films of the ma- 
terial act as spin-ladders and show field-induced super- 
conducitvity. While controversy surrounding the work of 
Schon et. al. puts the existence of field-induced super- 
conductivity in doubt, the fact remains that this material 
is very close to being quantum critical and could go into 
an ordered or disordered phase with small changes in pa- 
rameters. 

Here we consider two square-lattice layers of spins, 
with axes of the square-lattice pointing along x and y 
and a Heisenberg Hamiltonian 



a,i a.i i 

(1) 

Here, a takes values 1 and 2, the first sum runs over the 
nearest- neighbors along the x-axis, the second over the 
nearest neighbor along the y-axis and the third between 
the nearest-neighbors along the z-axis. The model 7i, in 
the limit of a = 1 (the isotropic bilayer), has been exten- 
sively studied using numerical and analytical methods— 
We shall call Ti. an anisotropic bilayer model. We set 
J — \ and study the phase diagram of this model at 
T = using Series Expansion Methods and Stochastic 
Series Expansion (SSE) Quantum Monte Carlo (QMC) 
simulations. We find that the phase boundary separat- 
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closely the behavior 

Jl = A^, (2) 

where A k, 2.53 is the critical point for isotropic case. We 
also calculate the uniform susceptibility and the internal 
energy of the model at finite temperatures, and study 
its excitation spectra in the spin-gap and Neel ordered 
phases. 

The plan of the paper is as follows: In section II, we 
discuss the series expansion method. In section III, we 
discuss the stochastic series expansion calculations. In 
section IV, we present the T — Q properties of the system. 
In section V, we present the uniform susceptibility and 
internal energy at finite temperatures and in section VI, 
we present our conclusions. 



II. SERIES EXPANSION METHODS 

We have carried out dimer expansions and Ising ex- 
pansions for this system at T = 0. The linked-cluster 
series expansion method has been previously reviewed in 
Ref. 0, and will not be repeated here. 



A. Dimer Expansions 

In the limit that the exchange coupling along the rung 
Ji_ is much larger than the couplings within the plane, 
the rungs interact only weakly with each other, and the 
dominant configuration in the ground state is the product 
state with the spin on each rung forming a spin singlet. 
We can construct dimer expansion in J/Jj_ by treating 
the last term in Eq. as the unperturbed Hamiltonian 
and the rest of terms as a perturbation. 

We have carried out the dimer expansions for the T = 
ground state energy per site, Eq/N, the antiferromag- 
netic susceptibility x, and the lowest lying triplet exci- 
tation spectrum A{kx, ky) (odd parity under interchange 
of the planes) up to order {J/J±y^ for fixed values of a. 
The resulting power series in J/ J±_ for the ground state 
energy per site Eq/N and the antiferromagnetic suscep- 
tibility X for a = 1/2 are presented in Table I. A table 
of series coefficients for the triplet excitation spectrum 
A(fca;, ky) would require an inordinate amount of space to 
reproduce in print; it is available from the authors upon 
request. Instead, we also present in Table I the series 
for the minimum energy gap m = A(7r,7r). The dimer 
expansion for isotropic case (a = 1) was carried out by 
Hida^^ up to order (J/Jj_)^ in 1992, and extended to or- 
der (J/J_l)^ by Gelfand^, and to order (J/J_l)^^ by one 
of the authorsSi. Here, the number of clusters involved 
is much more than the isotropic case since the system 
no longer has 90° rotation symmetry, and there are, in 
all, 38070 linked clusters of up to 12 sites involved in the 



B. Ising Expansions 

To construct an expansion about the Ising limit for this 
system, one has to introduce an anisotropy parameter x, 
and write the Hamiltonian for Heisenberg-Ising model as 

H = Ho + xV , (3) 

where 

a=l,2 j: a=l,2 i 

i a^i 
a=l,2 i 

+aJ ^ ''^{Sa,i^a,i+y + Sa,i^a,i+y) 
Q=l,2 i 

J2{Sl,Sl, + S^Sl^) - i ^ e„,,5^_, , (4) 

and ea,i = il on the two sublatticcs. The last term 
in both Ho and is a local staggered field term, which 
can be included to improve convergence. The limits x = 
and X = 1 correspond to the Ising model and the isotropic 
Heisenberg model, respectively. The operator Hq is taken 
as the unperturbed Hamiltonian, with the unperturbed 
ground state being the usual Neel state. The operator V 
is treated as a perturbation. It flips a pair of spins on 
neighboring sites. 

The Ising series have been calculated for the ground 
state energy per site, Eq/N, the staggered magnetization 
M, the uniform perpendicular susceptibility x±i ^nd the 
lowest lying triplet excitation spectrum A{kx, ky,kz) for 
several ratio of couplings and (simultaneously) for sev- 
eral values of t up to order x^'^ (the series for uniform 
perpendicular susceptibility x± is one order less). The 
series are available upon request. Here there are two 
branches of the spin- wave dispersion. From the series one 
can see that the symmetric (optical) excitation spectrum 
A.{kx,ky,0) is related to the antisymmetric (acoustic) ex- 
citation spectrum A(kx, ky, n) by 

A(fc^, fcj^,0) = A(7r - fcy,7r - fca;,7r) , (5) 

and so we only consider the antisymmetric excitation 
spectrum here. 

III. QUANTUM MONTE CARLO 
SIMULATIONS 

We have used the stochastic series expansion (SSE) 
method^SiS^ to study the ground state and finite temper- 
ature properties of the Heisenberg antiferromagnet on 
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FIG. 1: The ground state phase diagram in the parameter 
space of the in-plane anisotropy, a, and the inter-layer cou- 
pling, Jx. Both QMC and series expansion results are shown. 
The error bars are smaller than the size of the symbols. The 
curve is a square root fit to the critical coupling data. 



QMC technique based on importance sampling of the 
diagonal matrix elements of the density matrix e^'^^. 
Ground state properties are obtained by using sufficiently 
large values of (3. There are no approximations be- 
yond statistical errors. Using the "operator-loop" cluster 
update the autocorrelation time for the system sizes 
we consider here (up to « xlO'^ spins) is at most a few 
Monte Carlo sweeps even at the critical couplingi^ 



IV. T = PROPERTIES 
A. Results from Dimer Expansions 

With the dimer series for the antiferromagnetic sus- 
ceptibility X, and the minimum triplet gap to, one can 
determine the critical point {J/J±)c by constructing Dlog 
Pade approximanta^ to these series, and since the tran- 
sition should lie in the universality class of the classical 
d = 3 Heisenberg model (our unbiased analysis also sup- 
ports this), we expect that the critical index for x ^-nd 
TO should be approximately 1.40 and 0.71, respectively. 
The critical line obtained by the exponent-biased Dlog 
Pade approximantsSS are shown in Fig. ^ 

The spectra for some particular values of a and J/ Ji_ 
in the dimer phase are illustrated in Fig. [51 where the di- 
rect sum to the series at J/ J± = {J / Ji_)c is indeed consis- 
tent with the integrated differential approximant&^° that 
one can construct. In Fig. O we show the spectra for 
some particular values of a and J / J±_ along the critical 
line. 

To compute the critical spin- wave velocity w, one ex- 
pands the spectrum A in the vicinity of wave-vector 
(tt, tt) up to k^: 



and it is easy to prove^" that the critical spin-wave ve- 
locity along X and y direction are equal to {2CDxY/'^ 
and (2CDy)-^/^, respectively, at {J/J±)c- The series co- 
efficients for 2CDx and 2CDy in J/J± are listed in Ta- 
ble I. These series can be extrapolated to {J/J±)c by 
using the integrated differential approximants, and the 
results are shown in Fig. 0] one can see that Vx (vy) is 
increased (decreased) once the anisotropy is introduced. 
As a — > 0, the system approaches a spin-chain along the 
x-axis, where the spin-wave velocity is known to be tt J/2. 
Hence asymptotically, along the critical line, 



Vx/J± ~ (7r/2)(J/Ji) ^ 0.619/\^ 



(7) 



This asymptotic behavior is shown by a dashed line. One 
can see that the presence of other interactions further 
increases the spin-wave velocity along x. To understand 
the behavior of Vy, we turn to linear spin- wave theory. 
We obtain the antisymmetric excitation spectrum 



A{kx,ky) ^2SJ {l + a + J±/2jy 



cos 



(kx) + acos{ky) — Jj_/2J)' 



1/2 



(8) 



with this one can obtain spin-wave velocity along the y- 
direction 



Vy = 2S'J[a(l 



J±/J)] 



1/2 



(9) 



Note that Vx is obtained from Vy by exchanging the in- 
tralayer couplings along x and y-directions. The linear 
spin-wave results for both Vx and Vy are plotted as dot- 
ted lines in Fig. 0] We see that for small a, the linear 
spin-wave theory, which is known to be off by a factor of 
7r/2 for the spin-half chain, works quite well for Vy, the 
velocity along the weakly-coupled direction. 



B. Results from Ising Expansions 

With the Ising series for staggered magnetization M 
and the uniform perpendicular susceptibility x± , one can 
also determine the phase boundary by extrapolating the 
series for AI and x± to the isotropic point x — 1 using the 
same method as in Ref. 21, the results for a = 1/2 are 
shown in Fig. [S] We note that Af and x± first increase 
for small J±/ J, then decrease for a larger value of J±/ J, 
and vanish at about J_\_/ J — 1.8, which is consistent with 
the more accurate critical point determined by the dimer 
expansions. The reason for the initial increase is that 
for small J_l/ J the interlayer coupling enhances the anti- 
ferromagnetic long-range order as the system acquires a 
weak three dimensionality and quantum fluctuations are 
suppressed. 

The antisymmetric excitation spectra /S.{kx , ky , t:) for 
some particular values of J^/J and a = 1/2 are illus- 
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FIG. 2: Plot of the antisymmetric spin-triplet excita- 
tion spectrum A{kx, ky,TT)/ J± in dimer phase along high- 
symmetry cuts through the Brillouin zone for the system 
with couphng ratios a = 0.5 and J/J± = 0.1 (dotted hne), 
0.3(dashed hne), 0.556 (sohd line). The lines are the estimates 
by direct sum to the dimer series, and the points circles with 
error bar for the case of J/ J± — 0.556 only arc the estimates 
of the integrated differential approximants to the dimer series. 
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FIG. 4: The spin- wave velocity v/ J± along x and y directions 
versus a at the critical ratio {J/J±)c obtained from dimer 
series expansion. Also shown are the asymptotic results of 
Vx/J± = 0.619/\/a (dashed line) and the linear spin-wave 
results for Vx and Vy (dotted lines). 
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FIG. 5: The staggered magnetization M and the uniform 
perpendicular susceptibility x± vs J±/J for a = 0.5 as esti- 
mated by Ising expansions. 



FIG. 3: Plot of the antisymmetric spin-triplet excitation 
spectrum A{kx, ky,TT)/ J± along high-symmetry cuts through 
the Brillouin zone for the system along the critical line with 
a = 0.25 (dashed line), 0.5 (dotted line), l(solid line). The re- 
sults of the integrated differential approximants to the dimer 
series are shown. 



gapless at the (tt, 7r,7r) point. For Jj_/J = 1.75 (close 
to the critical point), the spectra are very similar to 
that obtained from dimer expansion. Thus, as one goes 
through the quantum phase transition, the spectra evolve 
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FIG. 6: Plot of the antisymmetric spin-triplet excitation 
spectrum A{kx,ky,-K) / J± in Neel order phase derived from 
the Ising expansions along high symmetry cuts through the 
Brillouin zone for the system with coupling ratios q = 0.5 and 
J_i/J = 0.5 (upper curve), 1.75 (lower curve), also shown as 
dotted lines are the results of linear spin-wave theory. 



C. Ground State Properties Obtained by SSE 

For the ground state properties, using the SSE method, 
we have studied lattices of the form 2 x L x L, with L 
up to 20. Periodic boundary conditions were applied in 
both the X- and y-directions. A series of values of a were 
chosen, and the critical J_l was determined for each of 
them, mapping out the ground state phase diagram in 
the a — Ji_ parameter space. An inverse temperature 
/3 = 16i was found to be sufficient for the calculated 
quantities to have converged to their ground state values. 

An efficient way to determine the critical coupling for 
the spin gap transition is by studying the finite size scal- 
ing of the ground state (T'=0) spin-stiffness. The spin 
stiffness can be definec&Si as the second derivative of 
the free energy with respect to a uniform twist (p. At 
T = 0, the free energy is the same as the internal energy, 
and the expression for the spin stiffness takes the form 



where E{(j)) is the internal energy per spin in the presence 
of a twist. The stiffness can be related to the fluctuations 
of the "winding number" in the simulation£^^>iS2iS2i^ and 
hence can be estimated directly without actually includ- 
ing a twist. Since the twist can be applied parallel to the 
X- or y-axis, there are two different spin stiffnesses, px 
and Py, in the anisotropic system considered here. 

Finite size scaling analysis dictates that at the critical 




FIG. 7: The ground state spin stiffness times the linear system 
size, L, as a function of the inter-layer coupling in the vicinity 
of the critical point for square lattices with L = 6, 8, ■ • ■ , 20. 
The (negative) slope increases with increasing L. The intra- 
layer anisotropy is a = 0.5. The upper (lower) panel shows 
the data for stiffness along the x-(y-) axis. Error bars are of 
the size of the symbols or smaller. The curves are quadratic 
fits to the data. 
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FIG. 8: The finite-size dependent critical point as obtained 
from the crossing of the stiffness curves for successive system 
sizes with linear dimensions L and L-\-2. The estimate for the 
critical coupling converges from above (below) for the Py{px) 
data, providing upper and lower bounds for the true critical 
coupling in the thermodynamic limit. 
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FIG. 9: The finite-size scaling of the staggered, Q=(7r,7r), 
static susceptibihty in the vicinity of the critical coupling for 
a=0.5. For lattices larger than L=10, the data close to the 
critical coupling is found to fit a power law behavior with an 
exponent 2 — r] {rj = 0.035) that is consistent with the 3D 
Heisenberg universality class. 



Size as 



p{L) ^ L 



d-2-2 



(11) 



where d{= 2) is the dimensionality and z is the dynami- 
cal critical exponent. The transition is expected to be in 
the universality class of the 3D Heisenberg model - hence 
z = \. It follows from the above relation that in a plot of 
Lpct, (a — X, y) versus J_l, the curves for different system 
size should intersect at J^. In practice, it is found that 
the crossing point shifts monotonically with increasing 
system size - the intersection of the curves for successive 
system sizes (linear dimensions L and L + 2) give a finite 
size dependent estimate of the critical coupling, J^L), 
that converges to the true critical coupling, J^, at large 
L. Interestingly, the convergence is from opposite direc- 
tions for px and py - J1{L) — > from above (below) 
for py{px)^ This yields upper and lower bounds for the 
true critical coupling, leading to an improved estimate for 
J^. The results for a=0.5 are shown in figEI The upper 
(lower) panel shows the data for px{py). The curves are 
found to cross in the neighborhood of J± « 1.8. A plot of 
the finite size dependent critical coupling obtained from 
the crossing of successive system sizes obtained from the 
plot is shown in figlHI As discussed above, J1{L) is seen 
to converge toward « 1.79 from above (below) for px{py). 
From the data, we estimate the true critical coupling in 
the limit of infinite system size to be — 1.79 ± 0.005 
for a = 0.5. 

A similar analysis of the spin stiffness data for different 
values of a gives us the critical value of the inter-planar 



a. The resulting phase diagram is shown in fign] 

Independent estimate of the critical coupling can be 
obtained from the finite size scaling of the staggered 
structure factor and the corresponding susceptibility. 
The full two-plane static spin susceptibility is defined as 

[SIM -SIM 



At the critical point, the staggered, q = (tTjTt), suscepti- 
bility for finite size systems scale with the system size as 
a simple power law behavior determined by the critical 
exponent rj: 



X{L, Jl) ^ L^-r, 



(13) 



On a plot of ln(x) versus ln(L), the data should fall on a 
straight line with slope 2 — rj. In the spin gapped phase, 
the staggered susceptibility should go to a constant at 
large L, whereas in the Neel phase with long range an- 
tiferromagnetic order, it should diverge faster than ex- 
ponentially with L. Fig|51 shows the plot of ln(x) versus 
ln(L) for fixed a{= 0.5) and five different values of J± 
in the vicinity of the critical coupling. For large values 
of L, the data for the critical coupling is indeed found 
to fall on a straight line. The slope of the line yields 
T] « 0.035, in close agreement with its value for the 3D 
Heisenberg universality class. The deviation from power 
law behavior in both the spin gapped and Neel phases 
are in agreement with the above discussion. Thus the 
transition is found to belong to the universality class of 
the 3D Heisenberg model. This is found to be true for 
any finite value of a. 



V. FINITE TEMPERATURE PROPERTIES 

A. Uniform Susceptibility 

We start the presentation of the finite temperature 
properties with the uniform magnetic susceptibility de- 
fined as 



I 
N 



T.[^b + SL][Sl,+Sl^] (14) 



where N is the size of the system. For finite temperature 
properties, we have carried out the simulations on lat- 
tices with rectangular geometry, Lx x Ly, Lx Ly. This 
is to reduce finite size effects. It was shown by Sandvik— 
for a system of coupled Heisenberg chains, the finite size 
effects depended monotonically for rectangular lattices 
whereas for square lattices the behavior was less well be- 
haved. While such effects are expected for ground state 
properties also, our limitation in terms of computational 
power has restricted us to the use of square lattices. For- 
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FIG. 10; The uniform susceptibility as a function of temper- 
ature near the critical coupling for a=0.5. At low tempera- 
tures, the susceptibility depends linearly on T, in agreement 
with the prediction from field theoretic calculations. Above 
the transition, in the spin gapped phase, Xu scales to zero at fi- 
nite T, while below the transition, in the Neel phase, Xu scales 
to a finite value at T=0. At the critical point, Xu extrapolates 
to zero at r=0. The data shown are for a rectangular lattice 
with iV= 128x32. 
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FIG. 11: The energy as a function of the temperature for the 
Neel and the spin gapped phases. The data for J_l=1.85 is 
found to deviate from a pure dependence at low tempera- 
tures. 



B. Internal Energy and Specific Heat 



and lower bounds for the critical coupling, leading to a 
reliable estimate of J^. Simulations at finite tempera- 
tures require considerably less computer power - hence 
we are able to study rectangular lattices. In particular, 
lattices with aspect ratio — 4Lj, have been considered. 

Results from the study of the non-linear a model pre- 
dict that at the critical coupling, the uniform susceptibil- 
ity, x(2^)i is linear in T (at low temperatures) with zero 
intercept. The region of linear T dependence gives an 
estimate of the quantum critical regime. The finite-size 
effects in the estimates of T > data decrease rapidly 
with increasing system size - the difference in the esti- 
mate of Xu between lattices with = 64 and 128 is of 
the order of the magnitude of the error bars up to the 
lowest temperatures studied. Henceforth, the data are 
presented for lattices with Lx—128. FigllOlshow the tem- 
perature dependence of the Xu for three different values of 
the inter-planar coupling close to the critical value. The 
data shown are for a rectangular lattice with Lx=128. It 
was found that the difference in the estimates of Xu at 
low temperatures For all the values, the uniform suscepti- 
bility is linear in T over the range of temperature shown. 
For Jl, the intercept is approximately zero, while it is 
positive (negative) for J± < (>) J^. The results are con- 
sistent with the estimate of obtained from the stiff- 
ness data. While the linear T behavior at the critical 
coupling is found to hold for all values of the in-plane 
spatial anisotropy, its range decreases with decreasing 

1 _r _ 



Instead of working directly with the specific heat, we 
have studied the internal energy, of which the specific 
heat is the temperature derivative. This is driven by the 
practical consideration that the specific heat data (that 
can be estimated directly within the framework of the 
SSE method) gets noisy at low temperatures, while the 
internal energy data is largely free from such noise at all 
temperatures considered. The temperature dependence 
of the specific heat is easily obtained from that of the 
internal energy. An estimate of the temperature depen- 
dence of the internal energy is obtained by assuming spin 
wave dispersion at low energies. In the Neel phase, the 
dispersion is e{k) — ck, where c is the spin wave velocity. 
In the presence of a spin gap (A), the dispersion takes 
the form £{k) — VA^ -|- c^fc^. For a 2D system, the inter- 
nal energy per particle at low temperatures should take 
the form 

For the gapless dispersion, this gives ^ = (const.) x T^, 
which implies that the specific heat, C„ ~ T^. For the 
spin gapped phase, the internal energy expression reduces 
to 

^ ^ (const.) X r 

For A^ ^ T, the leading behavior of the specific heat 
is once again ^ T^. For large gap and/or low tempera- 
tures, this will turn into an activated behavior, coming 

c _ J J J 1 -rj-i-- i^.j 1 
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However, this region has not been accessible to our Monte 
Carlo simulations. Fig ll II shows the internal energy as a 
function of the temperature for two parameter sets cor- 
responding to the gapless and spin gapped phases along 
with the fit. This suggests that unless one goes down 
to very low temperatures, near the critical line it would 
be difficult to tell from the specific heat data whether one 
is in an ordered or a spin-gap phase. 



cal Heisenberg model. The T = and finite temperature 
properties of the model are studied, especially around the 
critical line. We find that the spin-wave spectra evolve 
smoothly through the transition. The spin-wave veloc- 
ity becomes highly anisotropic as a —^ 0. We hope our 
work would stimulate further measurements of spin- wave 
spectra for materials such as CaCu203. 



VI. CONCLUSIONS 
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TABLE I: Series coefficients for dimer expansions of the ground-state energy per site, Eo/NJ±, the antiferromagnet suscepti- 
bility x±- the minimum triplet energy gap rn/J±, and the critical spin- wave velocity along x and y-direction (as described in 
the text). Coefficients of (J/Jx)" for a = 1/2 up to order n = 11 are listed. 



n 


£o/J± 




m/J± 


2CD^/Jl 


2CDy/Jl 





-3.750000000x10" 


1 


1.000000000 


1.000000000 




0.000000000 


0.000000000 




1 


0.000000000 




3.000000000 


-1.500000000 




1.000000000 


5.000000000x10" 


1 


2 


-2.343750000x10" 


1 


7.125000000 


1.250000000x10" 


1 


0.000000000 


0.000000000 




3 


-1.054687500x10" 


1 


1.546875000x10^ 


-4.687500000x10" 


1 


1.750000000 


5.468750000x10" 


1 


4 


-2.270507813x10" 


2 


3.132617187x10^ 


-5.078125000X 10" 


1 


1.218750000 


1.757812500x10" 


1 


5 


8.184814453x10" 


2 


6.153011068x10^ 


-1.647949219x10" 


1 


-4.101562500x10-1 


-3.376464844x10" 


1 


6 


9.038543701x10" 


2 


1.188306387x10^ 


2.856903076x10" 


1 


-1.158325195 


-1.558227539x10" 


1 


7 


1.045632362X 10" 


2 


2.276296512x10^ 


7.159423828x10" 


2 


-1.660129547 


-8.045387268x10" 


2 


8 


-1.300574541x10" 


1 


4.326861726x10^ 


-8.715667129x10" 


1 


6.815385818x10-1 


7.272934914X 10" 


1 


9 


-1.566342926x10" 


1 


8.165576512x10^ 


-2.380170342 




3.975312367 


8.960766308x10" 


1 


10 


-2.075107419x10" 


2 


1.529763646x10^ 


-1.890960221 




2.580596185 


-2.333164847x10" 


1 


11 


2.349058435X 10" 


1 


2.850279307x10^ 


-1.916425429x10" 


1 


-2.162085212 


-1.298405600 





